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An advection-diffusion-limited dissolution model of an object being eroded by a two-dimensional potential 
flow is presented. By taking advantage of the conformal invariance of the model, a numerical method is introduced 
that tracks the evolution of the object boundary in terms of a time-dependent Laurent series. Simulations of 
a variety of dissolving objects are shown, which shrink and then collapse to a single point in finite time. The 
simulations reveal a surprising exact relationship whereby the collapse point is the root of a non-analytic function 
given in terms of the flow velocity and the Laurent series coefficients describing the initial shape. This result is 
subsequently derived using residue calculus. The structure of the non-analytic function is examined for three 
different test cases, and a practical approach to determine the collapse point using a generalized Newton-Raphson 
root-finding algorithm is outlined. These examples also illustrate the possibility that the model breaks down in 
finite time prior to complete collapse, due to a topological singularity, as the dissolving boundary overlaps itself 
rather than breaking up into multiple domains (analogous to droplet pinch-off in fluid mechanics). In summary, 
the model raises fundamental mathematical questions about broken symmetries in finite-time singularities of 
both continuous and stochastic dynamical systems. 

I. INTRODUCTION 

Interfacial growth processes, such as alloy solidification [ 1- 
4], electrodeposition [5, 6], and crystal formation [7, 8], are re¬ 
sponsible for a wide variety of complex natural patterns [9, 10] 
that emerge due to instabilities in the underlying equations for 
interface motion [11, 12]. Often, continuum models of interfa¬ 
cial growth exhibit finite-time singularities, whereby features 
of the interface, such as curvature, diverge after finite time. The 
formation of such singularities is indicative of a breakdown 
in the separation of scales between the macroscopic variables 
in the continuum model and the microscopic variables that 
have been ignored [13, 14]. In fluid mechanics, the presence of 
finite-time singularities has been extensively studied, and can 
frequently provide physical insight [15], such as identifying 
the existence of universal scaling behavior in the pinch-off of 
a column of fluid [16, 17]. 

A particularly good example of interfacial growth is 
diffusion-limited aggregation (DLA) [18], where a solid cluster 
of particles is grown from a bath of diffusing particles start¬ 
ing from a single static seed particle. Additional particles are 
introduced one by one far away from the cluster, carry out 
random walks until they adhere to the cluster upon contact, 
causing it to grow. Since a random walker is more likely to 
first meet an extremity of the cluster than an interior region, 
the extremities grow preferentially, leading to complex fractal 
clusters in discrete computer simulations of the model [18-21]. 

Growth processes related to diffusion-limited aggregation 
have also been studied in the continuum limit, whereby the 
steady-state walker concentration satisfies Laplace’s equation 
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outside the cluster, is zero on the cluster boundary, and tends to 
a steady concentration far away from the cluster, and the cluster 
boundary grows continuously and deterministically with its 
velocity proportional to the normal gradient of the diffusing 
concentration field. This limit is mathematically equivalent 
to the classical Saffman-Taylor problem [22] of viscous fin¬ 
gering in a Hele-Shaw cell without surface tension [23, 24]. 
This model is conformally invariant, which simplifies the anal¬ 
ysis and allows it to be studied in detail in two dimensions 
using conformal mapping [25], as first formulated in 1945 by 
Polubarinova-Kochina [26, 27] and Galin [28] for applications 
to oil recovery and water filtration in porous media. Continuous 
diffusion-limited growth is notoriously unstable, and perturba¬ 
tions in an object’s boundary progressively sharpen, eventually 
leading to the formation of cusps in finite time [29-31]. In 
viscous fingering, these finite-time singularities are regular¬ 
ized by surface tension [24, 32, 33], which leads instead to tip 
branching instabilities and the formation of fractal fingering 
patterns [23]. 

Hybrid discrete-continuous models have also been devel¬ 
oped based on iterated conformal maps [34], which take full 
advantage of conformal invariance in two dimensions [25]. 
The growing cluster is defined by a chain of conformal maps 
that each add a small bump to the shape to represent the ag¬ 
gregation of a single particle, which opens new possibilities, 
such as growing non-random, fractal clusters [31]. By study¬ 
ing the statistics of the stochastic conformal map, it can also 
be shown that the average shape of random DLA clusters is 
similar [35], but not identical [36, 37], to the corresponding 
shape of continuous, deterministic diffusion-limited growth. 

Both the stochastic and continuous growth models have been 
extended to allow for any gradient-driven transport process in 
two dimensions [38], on flat or curved surfaces [39], which 
is made possible by a general conformal invariance princi- 
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pie [40]. In the canonical case of advection-diffusion-limited 
aggregation (ADLA) of random walkers in a fluid potential 
flow, asymptotic approximations of the flux profile have been 
studied [41], and the discrete and continuous cases have been 
compared [37]. These studies of ADLA provide the motivation 
for the present work. 

In this paper, we consider when the sign of growth is 
switched in the ADLA model, corresponding to dissolution or 
erosion: we start with a solid object, and then random walkers 
in a flow annihilate small parts of it on contact. Even in the 
absence of flow, this case has received much less investigation, 
since it usually leads to stable dynamics [21, 35, 42, 43] and 
thus many of the complex patterns due to growth instabilities 
are no longer manifest. However, this model opens up alter¬ 
native questions for study. In a previous paper [44], several 
different conformally invariant transport-limited dissolution 
models were introduced, including the erosion of corrugations 
on an infinite surface, and the expansion of a cavity due to 
dissolution. The paper also introduced the system of advection- 
diffusion-limited dissolution (ADLD), whereby an object is 
dissolved due to a concentration of random walkers in a fluid 
potential flow past the object. The object is represented by 
a time-dependent conformal map from the unit circle to the 
physical domain. By making use of previous asymptotic re¬ 
sults [41], an evolution equation for the conformal map is 
derived for the regime of intermediate Peclet number. The 
mathematical model is similar to the model of freezing and 
melting of dendrites in a hydrodynamic flow considered by 
Kornev and coworkers [45-47], where the random walker con¬ 
centration held in our model is replaced by a temperature held 
governed by an advection-diffusion equation. However, these 
works use a different mathematical formulation, employing the 
Schwarz function [48] to model the time-dependent boundary, 
and consider a different model for interfacial motion. 

The previous smdy of ADLD [44] was entirely analytical, 
and thus only considered the simple shapes of a circle and 
ellipse. Here, we investigate this model in more depth, and 
develop a numerical implementation that can simulate the dis¬ 
solution of arbitrarily shaped objects. By simulating arbitrary 
objects, we are able to investigate the model in substantially 
more mathematical detail, particularly in relation to the forma¬ 
tion of several different types of finite-time singularity. 


A. Physical applications of the model 

While the model that we consider is for a mathematically 
idealized two-dimensional case, it is worth considering what 
situations it could be applied to. For free gas or liquid flow, 
the model is unlikely to apply. Laminar flow at high Reynolds 
number Re can be approximated as a potential flow with vis¬ 
cous boundary layer of width 1 /v^l^, but our model assumes 
that the diffusion layer of width 1 / which is much wider 
than any viscous boundary layer for the regime of intermediate 
Peclet number that we consider. For melting, this would imply 
a small Prandtl number Pr, while for dissolution this would im¬ 
ply a small Schmidt number Sc = Pe/Re. However, for most 
liquids Sc ^ Pr ^ 1, since momenrnm diffuses much more 


rapidly than heat or mass. For gases. Sc ~ Pr ~ 1 since the 
same collisional mechanism governs mass, momenrnm, and 
heat transfer, but this still violates the model assumptions since 
the viscous and diffusion boundary layers would have similar 
size. 

One area where the model may apply is for dissolution 
or two-phase flow in porous media, which has relevance to 
water transport in soils or to flow in oil reservoirs, as in the 
seminal papers [26-28]. In this case, the fluid can be modeled 
using Darcy flow, and the object would represent a solidified 
region within the porous medium undergoing dissolution or 
melting (as considered by Kornev and coworkers [45-A7]). A 
similar model of nonlinear advection-diffusion in a potential 
flow is also applicable to viscous gravity currents in Hele- 
Shaw cells or porous-media flows, in which a heavier fluid 
spreads diffusively by gravity on an impermeable surface, as 
it is sheared by the flow of a lighter fluid above it [49]. The 
model could also have applications to a variety of different 
electrochemical corrosion processes [50-52], whenever a flow 
is imposed to modify the diffusive transport of active ionic 
species. 

The model may also be relevant to dissolution driven by 
electrokinetic phenomena, in the regime where double layers 
are thin in comparison to the object size. Consider a fixed 
object that is uniformly charged, in a fluid that is driven by 
a uniform electric held. In this case the fluid motion will be 
well-modeled by a potential flow where the electric potential 
is proportional to the fluid potential. The model could also 
apply to an object moving via electrophoresis. An object of 
constant surface charge will move at a constant speed in a 
uniform electric held, regardless of its shape [53], and hence 
constant flow at infinity could be fictitious flow of a stagnant 
fluid in a frame of reference of a particle moving at constant 
velocity by electrophoresis as it dissolves. 

The model we consider is in a substantially different regime 
than recent experiments and simulations of the erosion of clay 
bodies in high Reynolds number fluid flows [54], where the 
erosion rate of the surface is proportional to shear stress, and 
the dissolving bodies tend towards a self-similar cone-like 
shape [55]. 


B. Layout of the paper 

The paper proceeds as follows. In Section II, we present the 
theoretical background of the model, and derive an evolution 
equation for the shape of a dissolving body in terms of a time- 
dependent conformal map from the unit circle to the physical 
domain, described by a Laurent series. In Section III, we 
then derive a system of ordinary differential equations that 
govern how the Laurent series coefficients evolve with time. 
We numerically integrate this system using an eighth-order 
timestepping method [56, 57], which allows the dissolution 
process to be simulated very accurately, close to the limit of 
machine precision. 

Our initial numerical results for a variety of objects show 
that they completely dissolve in a finite duration with their 
boundaries becoming progressively smoother (Section IV). As 
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expected, the flow causes the objects to dissolve more quickly 
on the side facing upstream, although the details of the process 
are complicated, and affected by the precise manner that the 
fluid flows past the dissolving object. Of particular interest is 
the location of the collapse point, where the dissolving object 
finally vanishes. Due to the high accuracy of our simulations, 
we inferred an exact relationship between the collapse point 
Zc expressed as a complex number, the speed of the flow, and 
the initial Laurent coefficients. The relationship is surprising, 
whereby Zc is the root of a non-analytic function P, the terms 
of which involve complicated products of Laurent series terms. 
While some of these terms share similarities with binomial and 
multinomial expansions, they are distinctly different, and we 
are unaware of any other problem in conformal mapping or 
elsewhere where they occur. 

In Section V we make use of residue calculus to derive the 
general form of P, using the numerical results as a guide. The 
complicated products of terms in P arise from the residue of 
a contour integral where several Laurent series are multiplied 
together. In general, the function P has multiple roots, thus 
creating ambiguity about which root is the collapse point, and 
in Section VI we consider three different example objects that 
highlight the structure of P in more detail. To And the roots of 
P, we introduce a generalized Newton-Raphson iteration. As 
usual for Newton-Raphson iterations, plots of the root conver¬ 
gence in terms of the initial starting guess are fractal, but the 
non-analyticity of P creates some distinct morphological dif¬ 
ferences, and the plots illustrate the difficulties of determining 
the collapse point with mathematical certainty. 

The three examples also exhibit several types of finite-time 
singularity, which are both physically relevant and provide 
insight into the mathematical well-posedness of the model. 
Connections between these singularities and P are discussed. 
While the dissolution model that we consider is a simplified 
model with stable dynamics, it has a surprising amount of math¬ 
ematical structure, and our results raise a number of questions 
for further study. 


II. THEORETICAL BACKGROUND 

We make use of non-dimensionalized units, and consider 
an object in two dimensions with a time-dependent boundary 
S{t) as shown in Fig. 1(a). The object is immersed in an 
inviscid, irrotational fluid with velocity v(x,t), which can be 
written in terms of a potential 0(x,f) as v = V(j). The fluid is 
incompressible, so V • v = 0 and hence 

V2(j)=0. (1) 

At the boundary of the object the condition n v = n V0 = 0 
is used, where fi is an outward-pointing normal vector. Far 
away from the object the flow tends to a constant horizontal 
velocity so that v(x,t) —)■ (1,0) as |x| —>■ Equivalently, the 
potential satisfies ())(x,f) —>■ x as |x| —)■ 

The fluid transports a random walker concentration c(x,t) 
that satisfies the advection-diffusion equation 

PeVc • V0 = V^c, (2) 


where Pe is the Peclet number, a dimensionless quantity de¬ 
scribing the ratio of advection to diffusion. Far away from 
the object, the random walker concentration tends to unity, so 
that c(x,f) —>^ 1 as |x| —> o°. The random walkers are respon¬ 
sible for dissolving the object. At the boundary of the object, 
c(x,0) = 0. The normal velocity of the object boundary S(t) is 
given by 


a = —An-Vc, (3) 

where A is a dimensionless constant. Equations (1) and (2) 
together with the associated boundary conditions form a closed 
system for (()),c,5) that describe the dissolution dynamics, but 
they are difficult to solve directly. To proceed, we therefore 
treat the object as being in the complex z plane, where z = 
X + iy, and we introduce a time-dependent conformal map 
described by an analytic function z = g{w,t) that transforms 
the unit circle C into the object boundary S{t), as shown in 
Eigure 1(b). The most general form of the conformal map is 
the truncated Laurent series. 


8{w,t) 


N 


= a{t)w + Y, qnit)w ", 

n=0 


(4) 


where a{t) is taken to be a real function, and qn{t) are complex 
functions. Hereafter, we refer to q„ as the nth mode. Both 
Eqs. (1) and (2) are conformally invariant. The Laplacian is the 
standard example of a conformally invariant operator, and the 
advective term Vc • V()) is also conformally invariant [38, 40]. 

The boundary conditions in the w plane are different. Due 
to the scaling factor a(t) in Eq. (4), the boundary condition on 
the velocity potential becomes ^{w,t) —>■ aRe(w) as |w| —>■ 
We therefore introduce a rescaled potential = ^{w,t)/a 

that satisfies the original boundary condition ())(w,f) —Re(H’) 
as |w| —)• The rescaled system for c and satisfies Eqs. (1) 

& (2), but with a rescaled Peclet number Pe(f) = Pea(f). In 
addition, the normal growth in the w plane is 0^, — CT/|g'| to 
take into account the local volumetric scaling of the conformal 
map. 

Even in the w plane where the object is the unit circle, the 
concentration c cannot be determined analytically. However, 
asymptotic expansions have been studied in detail [41], and for 
Peclet numbers below 0.1, the approximation 


A/o(Pe)cP‘=“^e ^ 

<Jw ~-APe cos 0 


rPe 


h {t)e‘ 


J COS0 


^o(l) 


-dt 


(5) 


is uniformly accurate in 0 = argw. Taking the leading term of 
this approximation gives 


(Tw 


A(1 -fPe COS0) 
- 7 -log^ 


— APe cos 0, 


( 6 ) 


where 7 is Euler’s constant. 

To make progress, we now focus on the intermediate regime 
starting at small Peclet number and endingjirior to collapse, in 
which it is reasonable to assume that logPe is a constant [44]. 
By rescaling the time, we choose A = 7 -I- log ^ without loss 
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of generality. If the constant B = PeA is introduced, which we 
subsequently refer to as the flow strength, then Eq. (6) becomes 


Gw =—I+Ba{t)cos 9. (!) 


To transform this back into the physical domain, consider a 
point onthez(f) =g(w(f),f) on the boundary 5(f) of the object. 
Taking a time derivative gives z = g'w + g. Multiplying by wg' 
and taking the real part gives 

Re(wg'z) = Re(wg'g'w) +Re(wg'g). (8) 


Since the point in the w plane mapping to z(t) must lie on 
the unit circle, it follows that ww = 1 and hence Re(vi>w) = 0, 
so the first term on the right hand side of Eq. (8) vanishes. 
The motion of the point in the z plane is z = an where n is 
the normal vector written as a complex number. Taking into 
account rotation and scaling, the normal vector is given by 


g' 

l^'l kl 


(9) 


and hence the left hand side of Eq. (8) is 


Re(wg'z) = Re = Re(|g'|a) = (10) 

Combining Eqs. (7), (8), and (10) yields 


Eq. (11) becomes 


( 


N 

l+Ba cos 0 = Re 

ae 

- ^n{b„-ic„)e‘'‘^ 

V 


n=0 


n=0 

Eq. (14) is real, and can be expressed in terms of components 
cosnO and sinnO for n = 1,... ,77 + 1, plus a constant term. 
Equating both sides of the Eq. (14) in each component leads to 
2N + 3 coupled ordinary differential equations for the 2N + 3 
variables a(t), b„(t), and c„(f). Hence, other than for cases 
where these equations are degenerate, s will be uniquely de¬ 
termined in terms of s. Furthermore, since Eq. (14) does not 
feature any higher harmonic of sine and cosine, it follows that 
s exactly represents the time-evolution prescribed by Eq. (11): 
if a shape initially is described in terms of a Laurent series 
using terms up to q/i/, it will remain perfectly described by this 
Laurent series throughout the whole dissolution process. 

The details of equating each component of Eq. (14) are given 
in Appendix A. Equating the constant terms gives 

N 

ad-Y^n{bnbn+c„c„) = -l. (15) 

n=Q 

Equating the terms with factors cos(A-|- 1)0 and sin(A-|- 1)0 
gives 



Re{wg'g) =—1+Ba{t)cos 9, (11) 


abi\i = dNb]^, 


acN = oNcn, 


(16) 


which describes the dissolution process in terms of a time- 
dependent conformal map. For B = 0 it becomes the 
Polubarinova-Galin equation [26-28], which has been used 
in previous continuum DLA (viscous fingering) studies with¬ 
out advection [23, 29]. The incorporation of the Ba{t) cos 9 
term [44] represents the simplest extension to account for the 
general effect of advection [38] and is therefore a useful and 
interesting model to study in its own right. 


respectively. Equating the terms with a factor of sinn0 for 
n — gives 


0 = -d{n - l)c„_i +ac„-i 

N—n 

m=Q 

bmCm+n) 


(17) 


III. NUMERICAL METHOD AND IMPLEMENTATION Frnally, equating the terms wrth a factor of cosn0 for n 

1,... ,77 gives 


A. Discrete formulation of the governing equation 

We now make use of Eq. (11) to formulate a numerical 
solution technique. We represent the dissolving object via 
the time-dependent conformal map in Eq. (4) with a fixed 
value of 77 > 1. We write the qn{t) in component form as 
bn{t) + ic„{t), and describe the shape of the object by the real 
vector s(f) = {a,bo,bi,. .. ,bN,co,ci,. . .,cn), with a total of 
277 -f 3 components. Using the two expressions 

_ N 

Wg' =aw-Y, n{bn - (12) 

«=0 
N 

g = d+Y,ibn + ic„)w^", (13) 

n=0 


J3„ = -d{n-l)bn-i +abn-i 

N—n p 
m=0 

mdm-Vn) 


(18) 


where j3„ = Ba if n = 1, and )3„ = 0 otherwise. The combina¬ 
tion of Eqs. (15), (16), (17), and (18) can then be expressed as 
a linear system 


M(s)s = v(s) (19) 

where M and v are matrix and vector functions of s, respec¬ 
tively. By writing Eq. (19) as s = (s)v(s), the system can 

be integrated numerically. 
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(a) 



(b) 



FIG. 1. (a) The physical problem considered, where a two-dimensional object with time-dependent boundary 5(f) is dissolved a chemical 
concentration being transported by an incompressible potential flow, (b) A reference domain of the same physical problem but where the 
boundary is the unit circle C. A time-dependent conformal map z = g{w,t) describes the transformation between the two domains. 


B. Numerical implementation 

The simulations are carried out using double-precision float¬ 
ing point arithmetic, using LAPACK [58] to invert the lin¬ 
ear system in Eq. (19). To time-integrate the equation, the 
DOP853 integration routine described by Hairer et al. [56] 
is used. This routine uses the eighth-order, thirteen-step 
Dormand-Prince integration method that has the first-same- 
as-last (FSAL) property, requiring twelve function evaluations 
per timestep [59]. As described in more detail later, the com¬ 
ponents of s can sometimes vary rapidly, particularly close to 
the time of collapse. The DOP853 routine employs adaptive 
timestepping, which can retain accuracy in this situation. The 
routine estimates the local error [60] using a combination of 
fifth-order and third-order embedded numerical schemes. For 
all of the subsequent results, the timestep size Af is continually 
adjusted so that the absolute local error per timestep remains 
below a tolerance of 10^'^. If the estimated error of a timestep 
exceeds the tolerance, then the timestep is rejected and the 
integrator tries again with a reduced Af. 

There are three scenarios where the DOP853 integrator ter¬ 
minates early: (i) if a maximum number of timesteps is reached, 
(ii) if the equations are detected as stiff [57], or (iii) if the 
timestep Af required achieve the desired local error becomes 
too small. In the following results, we have only observed 
the third scenario. This occurs when Af becomes smaller than 
lOMrf, where = 2.3 x 10^^® is an estimate of the smallest 
number satisfying 1.0 -F Mr > 1 0 in double-precision floating 
point arithmetic. In certain cases, such as the examples of Sub¬ 
secs. VIA and VIC, the third scenario signifies a breakdown 
of the physical problem due to the formation of a cusp. How¬ 
ever, the third scenario also occurs in normal cases close to the 
time of collapse f^ due to a(f) varying rapidly. If the DOP853 
integrator terminates within 10^Mr of tc then we manually ad¬ 
vance to tc using timesteps of lOMr or less. While this may no 
longer achieve the required level of local error, we find that it 
provides several additional digits of accuracy in the collapse 
point location, which is useful in some of the later analysis. 

In some of the subsequent results, we must evaluate s at time 
points spaced at fixed intervals, which may not precisely coin¬ 


cide with the time points that are selected during the adaptive 
time-integration, which are usually unevenly spaced. To solve 
this we use the dense output formulae described by Hairer et 
al. [56]. By performing three additional integration steps, the 
solution can be approximated as a seventh-order polynomial 
over the interval of a timestep, allowing s to be evaluated at 
any specific time point. For computational efficiency, these 
three additional steps are only done when one or more output 
time points overlaps with the current timestep interval. 

The simulations are implemented in C-H-i-, and the code 
required to perform all of the subsequent analysis is provided 
as Supplementary Information. For all of the results presented 
here, the computation time required to simulate the dissolution 
process is negligible, taking less than 0.25 s on a Mac Pro (Fate 
2013) with an 8-core 3 GHz Intel Xeon E5 processor. 


IV. RESULTS 


A. Analytic results for the area and highest mode amplitude 


Before presenting results of the numerical method. It is use¬ 
ful to establish some basic features of the equations presented 
in the previous section. The area of the object is given by 

A{t) = JJ^dz (20) 

where Q. is the region enclosed by S{t). Using Green’s identity 
in complex form, 

A{t) = -l-.i zdz=\-. (g{w)g’{w)dw. (21) 

2i Js(t) 2i Jc 

Since ww = 1 on the unit circle, the integrand can be converted 
into an analytic function. 




+ qnw" 
n=0 


dw 

( 22 ) 
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and applying residue calculus gives 

A{t) = 7i:(a^- = 7Z (- Y, n{bl + clU , 

\ «=0 / \ «=0 / 

(23) 

describing the area as a function of the current mode ampli¬ 
tudes. Furthermore, time-integrating Eq. (15) gives 

N 

a^-Y>^{bl+cf,)=C-2t, (24) 

«=o 

where C is a constant, and hence 

A{t)=Ao-27:t, (25) 

where Aq is the initial area of the object. The area of the object 
therefore decreases at a constant rate, independent of the flow 
parameter B, with the time to collapse given by 

The modes in Eq. (16) also have first integrals, 

bN = ka^, CN = la^ (27) 

for some constants k and /. The highest mode amplitudes are 
therefore only dependent on the conformal radius a. Due to the 
couplings in Eqs. (17) and (18), similar results for the lower 
modes do not exist. 


B. Initial numerical results 

Figure 2 shows the dissolution process for six objects calcu¬ 
lated using the numerical code, where for all cases a{Q) = 1 
and B = 0.7. Figure 2(a) shows the dissolution process for 
a circle. Throughout the process, the circle retains its shape 
although its center progressively moves rightward due to the 
effect of the flow, which preferentially dissolves the side of 
the circle facing upstream. A similar behavior is visible in 
Fig. 2(b) for an ellipse, which keeps its shape through the 
dissolution process, while the ellipse center moves up and 
right. The results for Figs. 2(a) and 2(b) match those that were 
previously studied analytically [44]. For the case of no flow 
where B — 0 where the system is mathematically equivalent to 
(time-reversed) Laplacian growth, and bubble contraction in a 
porous medium, it is known that the ellipse is the most generic 
self-similar shape [61-63]. 

Figure 2(c) shows the dissolution process for a triangular¬ 
shaped object given by setting q 2 = —0.35 initially. In general, 
the mode is responsible for an (n + l)-fold perturbation of 
the boundary. If all of the q„ are initially real, the object is sym¬ 
metric about the x axis, and will remain symmetric throughout 
the dissolution process. For the case shown, the point of the 
triangle that faces upstream is more rapidly dissolved than the 
other two. Unlike the previous two examples that retain their 
shape during dissolution, the triangle becomes progressively 


more rounded at later times. Figure 2(d) shows the dissolution 
process when the previous object is rotated by 90°, which is 
achieved by setting q 2 = 0.35/. This object is initially symmet¬ 
ric about the y axis, but the flow causes this symmetry to be 
lost as time passes. The collapse point is slightly up and right 
from the origin. 

Figure 2(e) shows the dissolution process for the case when 
^15 = 0.05/ initially, which creates a 16-fold perturbation in the 
boundary. After 20% of the object has dissolved, this pertur¬ 
bation is almost completely removed, with the object’s shape 
approaching that of a circle. This is expected from Eq. (27), 
which shows that the highest mode will be proportional to 
and hence decay more rapidly for larger A. If several modes are 
initially non-zero as in Fig. 2(f) an irregular shape is formed, 
which behaves like a combination of the previous examples, 
with sharp features in the boundary being rapidly removed. 

We now examine the evolution of the modes and look in 
detail at the effect of the flow strength B. We make use of 
the specific example of a diamond shape given by a = 1 and 
q-} = 0.25 initially. Figure 3(a) shows the dissolution process 
for the case of zero flow when B = 0. Similar to Figs. 2(c) 
and 2 (d) the object becomes progressively more circular, but 
without the presence of flow it retains symmetry in the x axis, 
y axis, and the line x = y. Figure 3(b) shows the time-evolution 
of the modes throughout the dissolution process. The modes 
qo, qi, which were zero initially, remain zero throughout 
the dissolution process—this is expected since any non-zero 
contribution from these modes would break at least one of 
the symmetries seen in Fig. 3(a). The dissolution process is 
therefore described entirely in terms of a and ^ 3 , and could 
therefore be determined analytically using Eqs. (24) and (27), 
as considered in previous work [29, 44, 64]. Since q^ remains 
at zero, the collapse point is at the origin. 

Eigure 3(c) shows the dissolution of the diamond when the 
flow parameter is B = 0.7. As in the previous examples of 
Eig. 2, the diamond dissolves more rapidly on the side facing 
upstream, and the collapse point is slightly downstream. The 
time evolution of the modes (Eig. 3(d)) is significantly altered 
in this case, with all three components qQ, q\, and q 2 becoming 
non-zero during the dissolution process, due to the mixing 
between modes via the advection term in Eq. (11). The effects 
of these three modes, such as the translation of the object 
center, and the loss of symmetry about the y axis, are clearly 
visible in Fig. 3(c). The qi and q 2 modes decay to zero at the 
point of collapse, while the q^, mode remains positive. The 
value of q'o at f = U gives the collapse point position. 

Figures 3(b) and 3(d) also indicate the adaptive integration 
timesteps chosen by the DOP853 integration routine. In the 
middle of the dissolution process, at f « 0 . 2 , the routine is 
able to take timesteps up to approximately 0.02 while retaining 
the desired level of local error of 10^*^. However, close to 
t = tc, many more timesteps are needed to resolve the rapid 
change in a. For the example shown in Fig. 3(d), a total of 
329 integration timesteps are evaluated. During the DOP853 
integration routine, 194 steps are accepted, and 128 are rejected 
due to the local error estimate exceeding the given tolerance. 
Six additional small steps are required to reach the collapse 
time tc- 
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FIG. 2. Sample dissolution processes for six objects, all starting with a = 1 and using B = 0.7. The six objects and initial non-zero modes 
are (a) a circle, (b) an ellipse with qi = 0.3 + 0.2i, (c) a triangle with q 2 = —0.35, (d) a triangle with q 2 = 0.35i, (e) a corrugated circle with 
^15 = 0.05i, and (f) an irregular object with qi = —0.28 + 0.2; and q^ = 0.1. The white lines show the flow streamlines around the initial shape. 
The colored regions shown the shapes of the object at successive times as it dissolves, where each progressive region represents the dissolution 
of 20% of the object’s initial area. The black circles indicate the final points of collapse. 
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FIG. 3. (a) Dissolution process of a diamond shape initially described by non-zero modes a = 1 and q-^ = 0.25 for the case of zero flow, B = 0. 
The colored regions shown the shapes of the object at successive times as it dissolves, where each progressive region represents the dissolution 
of 20% of the object’s area, with the black circle indicating the point of collapse, (b) Time-evolution of the modes describing the diamond 
during dissolution, where the small circles on each curve show the integration timesteps using the adaptive DOP853 integration scheme, (c) 
Dissolution process of the diamond when the flow is B = 0.7. (d) Time-evolution of the modes describing the diamond during dissolution with 
flow, with the small circles on each curve showing the integration timesteps. 


C. Inferring analytic formulae for the collapse point 

Figures 2 and 3 show that the collapse point Zc =Xc + iyc of 
the dissolution process is dependent on both the flow strength 
B and the initial shape of the body as described by its Laurent 
coefficients. Since there are no other quantities in the problem, 
Zc must be given in terms of B and the Laurent coefficients 
only. The precise form of this dependence is not obvious, as 
the collapse point is given as the component q(, of the nonlinear 
differential equation system, evaluated at the time of collapse 
L- 

In this section, we infer the exact form of this relationship 
by exploiting the very high accuracy of the simulations, which 
allow the collapse point to be calculated to at least twelve deci¬ 
mal places. To simplify the analysis, we set a = 1 throughout 
this section. To begin, we restrict to the case when the Lau¬ 
rent coefficients are given purely in terms of real components 
bj. As discussed in the previous section, the components will 


remain real throughout the simulation, and the object will be 
symmetric about the x axis. Hence the collapse point Zc will be 
real, and determined entirely in terms of the horizontal position 

Xt. 

Figure 4 shows an example of inferring an analytic rela¬ 
tionship, for the case of an ellipse where the only non-zero 
Laurent coefficient is hj. Figures 4(a-d) show four dissolution 
process are shown for when b\ = 0.3 and the flow strength is 
B = 0,0.3,0.6,0.9, respectively. The collapse point Xc moves 
progressively right as B is increased. In Fig. 4(e), the nu¬ 
merically computed Xc is plotted as a function of B, for three 
different values of bi of 0.0, 0.3, and 0.6. The plot demon¬ 
strates that Xc is linear in B, with the slope depending on bi. 
The numerical data matches the relationship 


B{l+bi) 


(28) 


with the sum of square residuals being 1.3 x 10 1.3 x 
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FIG. 4. An example procedure to infer the analytic formula for position of collapse point in terms of the initial modes, (a-d) Dissolution 
processes of an ellipse given by a = i,q\ = 0.3 for flow strengths of B = 0,0.3,0.6,0.9, respectively, using the same visual representation as 
described in Fig. 2. (e) Plot of the numerically computed horizontal collapse point position Xc = Rezc as a function of B for three different 
values of i>i, which match linear relationships to numerical precision, suggesting an analytical relationship. 


and4.5 X lO^^'^forBi = 0,0.3,0.6, respectively. These 
small residuals, which are of a similar size to the expected nu¬ 
merical error, strongly suggest that this is an exact relationship 
for the original mathematical problem. 

One can extend this analysis to the case where the only non¬ 
zero Laurent coefficient is b„, and determine that Xc satisfies 
the polynomial relationship 

B Xc-bnX- 
2 l-nbl' 

For n = 1, this is consistent with Eq. (28) for the ellipse, 
although it also reveals more structure, and by substituting 
Eq. (26) the relationship simplifies to 

Btc=Xc-bnX"^. (30) 

To proceed, we now consider if there are two non-zero Laurent 
coefficients. The simplest case would be bo and b„ being non¬ 
zero, for « > 1. Since bo corresponds to a translation, the 
relationship is immediately given by 

Btc = {xc - bo) - bn{xc - boY, (31) 

without the need for simulation. Expanding the second term 
yields 

Btc = -bnx"^+nbnbox^^^ - '^^^^b„boX^^^ 

+ «(n-l)(«-2) ^^^3^-3 - + (^^ - ^o) (32) 

where the coefficients on the powers of Xc follow Pascal’s 
triangle. 


The next case to consider is when bi and b„ are non-zero, 
for n > 2. Unlike the previous case this cannot be immediately 
derived, and must be Inferred through fitting to simulation. Ta¬ 
ble I shows the derived results for the cases of n = 2,3,..., 10 
where a surprising pattern emerges. We see polynomials that 
bear some resemblance to a binomial expansion, although in 
contrast to Eq. (32), only every second power of Xc is present. 
Furthermore, the coefficients in front of the terms are integer, 
but of a more complicated form than Pascal’s triangle. Unlike 
the previous case, the more complicated form of these polyno¬ 
mials precludes rewriting them in a succinct form like Eq. (31). 
The pattern continues for the case when b 2 and b„ are non-zero, 
for n > 3. As shown in Table II, only every third power of Xc is 
present. The integer coefficients follow a natural progression 
from those in Table 1. 

In Tables I and II, we observe that each pair of non-zero 
Laurent coefficients leads to a combination of additional terms 
appearing in the collapse point polynomial. Building on these 
results, we inferred and numerically tested the formula 

Btc = —b4x‘l—xl +x^{4b4bi — B 2 ) 

+Xc{l —b\ + 3 B 3 B 1 -I- 4 B 4 B 2 ) 
4 - 2 ^ 2 ^ 1 + 3 ^ 3 / 724 - 4 ^ 4^3 — 2 ^ 4 / 7 ] (33) 

for the case of all four coefficients b\, b 2 , /’ 3 , and b 4 being 
non-zero. In Eq. (33) all terms involve powers of two different 
bn, but for higher non-zero Eaurent coefficients, terms with 
three or more different b„ arise. If bi, b 2 , and bs are non-zero, 
then we find that 

Btc = -bsx^c+^bsbix^c + ( 5 / 75/72 - b 2 )x^ 

+ {l-bi- 5 / 75 / 7 ?) + 2/72/7] - 4 / 75 / 72 / 7 ], (34) 
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n 

2 

3 

4 

5 

6 

7 

8 

9 

10 


Q{xc) _ 

-b 2 xl + 2 b 2 b\ + {l-b\)xc 

—b$xl +3b2biXc + (1 —bi)xc 

—b 4 x‘^ + 4-b4bix^ — 2 b 4 b\ -\-{l—b\)xc 

-b^xl + 5bsb\xl - 5bsb\)xc + {I - bi)xc 

—b^^ + 6b(,bix‘l — 9b(jb\x^ + 2b(,b\ + {\—b\)xc 

—bjxl + lbibix^^ — lAbjblx^^ + lb-ib\xc + (1 —b\)xc 

—b^x^ + ^b%b\x^ — 20b^blx^ + I6b^b\2^ — 2b^b\ + {l—b\)xc 

—bgx^ + 9bgb\x]. — 21bgb\xl + 30bgb\xl + 9bgb\xc + {I — bi)xc 

—bioxjP + lQb\Qb\x^^ — 35bu)b\x^^ A-5Qb\ab\x‘^ — 25b\Qb\x^ -\-2b\Qb\ + (1 — ^i) 

-b„x"+nbnb\x’l^^ - b„bjx"^‘^ + b„b\x''^^^ - ... + {\-b\)xc 


X 


C 


TABLE I. Examples of the analytic relationship Btc = Q{xc) for the horizontal collapse point position Xc that were inferred numerically using 
the high-precision calculations, for the case of an object described by two real non-zero Laurent coefficients bi and b„. The integer coefficients 
colored in blue, green, and red follow patterns. The final line of the table shows an inferred general formula. 


n 

3 

4 

5 

6 
1 
8 

9 

10 
11 
12 


Q{xc 


-b^xl 4- 3b3b2 + {xc - b2xl) 

-b4X^^ 4- 4b4b2Xc 4- {xc - b2xj) 

-b^xl -I- 5bsb2x] 4- {xc - b2xl) 

-b(,x^ 4- 6b(,b2xl - 3b(,b2 + {xc - b2xl) 

—bjxl -|- 7 bjb2x‘^ — Ibjb^Xc 4- {xc — b2X^) 

-b%xl 4- 8bsb2xl - 12^8^2-^? + i^c - b2X^) 

—bgx'^ 4- 9bgb2x^ — iWgb\x^^ 4- 3bgb\ + (x^ — b2^) 

-biQX^^ + \Qb\ 0 b 2 xl - 25feio^2-*c + 10/tiofe2-^<; -I- (xc - b2xl) 
—bi\x\^ 4- llb\\b 2 X^^ — 33^11^2-^^ 4-22b\\b\x^ 4- (x^ — ^2-*?) 
-bi2X^^ + 12fei2^2-*c “42^12^2-^® +40^12^2-^e “ 3/712^2 + (-^e 
—b„x" 4-nbnb2X^^^ — ”^”2 bnb^x"^'^ - 


-bixl) 

■ ■ ■ + (Xe 


bixl) 


TABLE 11. Examples of the analytic relationship Btc = Q{xc) for the horizontal collapse point position Xc that were inferred numerically using 
the high-precision calculations, for the case of an object described by two real non-zero Laurent coefficients ^2 tmd b„. The integer coefficients 
colored in blue, green, and red follow patterns. The final line of the table shows an inferred general formula. 


where the last term on the right hand side is a product of all 
three non-zero Laurent coefficients. 

The final generalization that we consider is when the Laurent 
coefficients are complex. The fitting procedure described in 
Fig. 4 becomes more complicated in this case, since both the 
horizontal position Xc and the vertical position jc of the collapse 
point will vary. For the case of the Laurent coefficients q\ and 
q 4 being non-zero and complex, we inferred the formula 


which is equivalent to the formula for an ellipse = f (1 + ) 

that was derived in previous work [44]. 


V. DERIVATION OF THE COLLAPSE POINT FORMULAE 


Btc =Zc- qiZc - q4Zc+4q4q\zl - 2q4q\, (35) 


which is a generalization of the formula for n = 4 in Table 1. 
The generalization to complex coefficients introduces conju¬ 
gates on some terms, and the right hand side is not an analytic 
function of Zc due to the first term featuring Zc- If q '4 = 0 also, 
then Eq. (35) simplifies to 


B _ Zc - q\Zc 

2 \-q\q\' 


(36) 


The previous section revealed a surprising and complicated 
connection between the collapse point, initial shape of the 
object, and the flow strength. Using these numerical results as 
a guide, we now analytically derive this connection. While the 
formulae in Tables I and II are complicated, it is reasonable 
to imagine that the specific coefficients could occur as the 
residue from a contour integral, perhaps involving the product 
of several Laurent series, and thus our first step is to consider a 
general integral quantity and determine its behavior during the 
dissolution process. 
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A. Time-evolution of an integral quantity 


This can be written as a contour integral as 


Consider the expression 

/(f) = I F{z)dz 
Js{t) 


(37) 


HI r^TT 

— = / / F'{g{e'^)){2-Bae'^-Bae-'^)dd 

at Jo 

f F'{g{w)){2w — Ba —Baw^)dw 
fc 


(46) 


where z = giw), S{t) is the shape of the object, and F is an 
arbitrary analytic function. This can be written as 

7(0 = j^F{g{w))g'{w)dw, (38) 

where C is the unit circle, and w = e'®. Since vi>w = 1 on the 
unit circle, this can be converted into the integral of an analytic 
function. 


yielding the second integral expression that will be used later. 
As a check, it is useful to consider when F'(z) = 1, in which 
case the integral matches the one in Eq. (21) for the area of the 
object. Then 


dl 

dt 


(2w — Ba — Baw^)dw 




2%i{2) =AKi 


(47) 


and 


= ( 39 ) 

which can be written as 

I{t)= £F{g{w))-^{g{^))dw (40) 

and hence integration by parts can be used to obtain 


7(f) = - qnw"^ dw 

= —27ti ^ n\qn'\^ = —2niA{t) (48) 

This gives A(f) = —2, which agrees with Eq. (24). 


t{t) = -£F\g{w))g\w)g{l,)dw. (41) 


B. A specific integral quantity 


This is the first of two expressions that will be used later. To 
obtain a second integral expression, consider taking the time 
derivative, which gives 

- j^F'{g{w)){g'{w)g{^)+g'{w)g{^))dw. (42) 

Applying integration by parts to the first integral, transfers the 
derivative onto the g{w)g{\ jw) terms. Note that 

=8M8{i)-^8{i) (43) 

and since the first term of this expression will cancel with one 
of the terms in second integral of Eq. (42), it follows that 

^ = (^^8'{i)+8'Mg{^,)^dw. (44) 

By substituting w = e'® and making use of Eq. (11), 

flj rln , _ 

-=-y^ F\g(^F^))i(g[<A^)g'[<d^)e-'^ 

= -2ij\\g{e‘^))Re (7(?^g(e‘®)) dO 

rlTt 

= 2i F\g{F^)){l-Bacos9)d9. (45) 


A interesting candidate for the function F' is 

F\z) = — (49) 

Z-Zc 

where Zc is the collapse point. This function is particularly 
special, since as the object is dissolving, the integrals given 
in Eqs. (41) and (46) will always be finite, as the integration 
contour will never pass over the singularity. Even though the 
function F{z) = log(z — Zc) is multivalued, the two integral ex¬ 
pressions that will be used in the following derivation, Eqs. (41) 
and (46), are related to each other through a derivation only 
involving F', and thus it is not necessary to consider branch 
cuts that would be needed to integrate F. Equation (46) will 
give 


dl f {2w — Ba — BaMp‘)dw 
dt fcw^{aw-Zc + 'L^=oqnW-'‘) 
jT {2w — Ba — Baw^)dw 

Since ^(w,f) is a conformal map that takes the region |w| > 1 
to the region outside the object, there can be no solutions to 
g{w,t) = Zc for |w| > 1 and thus the above integrand will have 
no poles for |w| > 1. Hence the integration contour can be 
deformed outwards and evaluated in terms of the residue at 
infinity, which is given by the coefficient of the term, 
namely —Baja = —B. Hence 

dl 

-^ = -2niB (51) 

dt 
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and therefore /(f) = D — 2niBt for some constant D. To de¬ 
termine the constant, we consider the limit as f —> tc, where 
a —> 0, qo —> Zc, and q„ —)• 0 for all n > 0. Then Eq. (41) shows 
that 


/(f) ^ 


f a{^+qo)dw 
Jc {aw-Zc+qo) 


i ‘^dw = -Inizc 
Jc aw 


(52) 


and therefore /(f) = 2Ki{B{tc — t) —Zc)- Returning to Eq. (41), 
and examining the integral at the initial time. 


m = -jc 


{a - Ilo ^ + lOLo 

aw-Zc +1,^=0 qnW-'‘ 


(53) 


which for large w can be expanded as 


/(O) = - 


{a - Lto^nnw-("+^^) dw 

iaW-Zc+Ln=0^nW-'') 

{a - Lto^nnw-('‘+'^) dw 


c aw 




- + Y, 

n=0 



yto 


(54) 


This integral will give the desired relationship between Zc and B. By expanding out the three power series, and looking at terms of 
the form that will give a residue at inhnity, the integral will simplify to a polynomial in Zc- For example, consider the case of 
only q\ and q^ being non-zero. In that case, for w large, and neglecting terms smaller than 


/(O) = - j (- +^iw-f 

Jc Vw aw^ aw^ J Vw 


\<r=0 “ 


_ 3 , _ q4q\w 

q4W +qi - 

C \ a 

; 7.,4 


i + l(“ 

a Vw 


1 ( g4Zc 5q4q\zl ^ q^q] ^ q\Zc 

a 


w ^ a4 

^„-f94Zc 4q4qizl 2q4q\ 

= -2%l — -3--f-^ 

V 


w^ 

q4q\zl 


(Zc 

qi 

?4 Y 

V w 

w^ 

/ 

?4\ 

1 

(Zc 

/ 


V w 


dw 


q4q\\ 
a2 ) 


w^ 
dw 


w 


,5 ) 


-f... \dw 


a J 


(55) 


By using /(O) = 2ni{Btc — Zc) it follows that 

a'^Btc = a^Zc — q4Zc +4aq4qizl — 2a^q4q\ — a^qiZc- (56) 
If a = 1 then 

Btc =Zc- q4zt +4q4qizl - 2q4q\ - qiZc, (57) 
which agrees with Eq. (35) that was found numerically. 

VI. THREE EXAMPLES OE THE COLLAPSE POINT 
EQUATION 


expansion procedure described in the previous section. The 
polynomial can then be analytically extended to give an ex¬ 
pression for P{zc) at points outside the object also. However, 
at points outside the object, the analytic extension will not 
match the value of the integral, since the enclosed residues will 
be different. Equation (58) is complicated: it is not analytic 
due to the presence of Zc, and in general it will contain higher 
powers of Zc, so it is likely to have multiple solutions. To use 
this equation as a predictive tool, it is useful to understand the 
typical structure of P and know how to select the correct root. 
We now consider three examples that explore the structure of 
P in relation to the object shape. 


For a general case, the collapse point Zc satishes the equation 

0 = P{zc) =Zc-Btc + :^.(f log(z - zc)dz, (58) 

27tl Js{ 0 ) 

defined at points within the object, where the integral in this 
equation is evaluated as a polynomial in Zc following the series 


A. First example: an irregular pentagonal shape 

Consider an example based on Eq. (57), where the function 
P can be written as 

PiZc) =Zc-Btc-q4zt+4q4q\Zc-2q4qi-qiZc, (59) 
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FIG. 5. (a) The thick black line shows an first example object, where a = 1 and the only non-zero Laurent coefficients are 
qA = TO + Ths colors show the argument of the function P[z) for B = ^, whose roots represent candidates for the collapse point of the 
object as it dissolves. The dashed lines are contours of |F’(z)| at the values of " for u e N. (b) A zoomed-in region showing forward and 
backward time-evolution of the object boundary at intervals of The unique negative-sense root of P is shown by a circle, and one of the 
positive-sense roots is shown by a triangle. The four other positive-sense roots are outside the region that is plotted. 


where = 1 — |?ip —4|q'4p. Fig. 5 shows a plot of the mod¬ 
ulus and argument of this function for the case of a = 1, 
q\ = ^ qA = ^ and B = ^. The shape of the 
object is also shown. There are five roots that lie outside the 
object. There is one root inside the object, which must be the 
collapse point. Furthermore, the argument in the neighborhood 
of the interior root rotates in the negative (anti-analytic) sense, 
whereas the argument near each exterior root rotates in the pos¬ 
itive (analytic) sense. By considering the Taylor series of P at 
a given root, one can mathematically determine whether a root 
is positive-sense or negative-sense by whether — \Pi^ is 
positive or negative, respectively. From Eq. (52), the collapse 
point must be given by a negative-sense root, and hence for 
this example there is an unambiguous choice, of the single 
negative-sense root within the object. 


M(s) becomes singular, and the DOP853 integrator terminates 
because the timestep required to keep the local error below 
the tolerance is smaller than what can be resolved with double 
precision. While the positive-sense root appears connected to 
the development of the cusp, it is not located exactly at the 
cusp, and thus it is not clear what, if any, its precise physical 
significance is. 

A practical way to determine the root positions is to make 
use of a Newton-Raphson iteration, generalized to take into 
account that P also depends on the conjugate of Zc- An ap¬ 
propriate Newton-Raphson iteration can be constructed by 
viewing P as a function of two variables Zc and Zc, and consid¬ 
ering the two-function system of P and P. For a guess of the 
form the vector generalization of the Newton-Raphson 
method to give an improved guess is then 


It is interesting to consider whether the other roots have 
physical significance. Figure 5(b) shows a zoomed-in region 
of the dissolution process for this example, confirming that 
the interior negative-sense root visible in Fig. 5(a) is indeed 
the collapse point. The figure also shows a nearby positive- 
sense root. If the system is time-integrated backward, then 
the boundary of the object sharpens toward the root. This 
leads to a cusp singularity in a finite time t = —0.06133, which 
appears similar to cusp development in related systems [24, 
29, 64]. The cusp formation occurs when a branch point of g 
reaches the unit circle. As the cusp is approached, the matrix 


f Pz P-z 
\Pz p-z 


(n+1) (n) 

Zc — Zc 


Jn+l) Jn) 
Zc Zc 



(60) 


which leads to the two equations 

-Z^c'’)+Pz{z^c^'^-Z^c^) = -P, (61) 

) + Piiz'f^^^ - ) = -p. (62) 


Substituting Eq. (62) into Eq. (61) to eliminate 
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FIG. 6. Plot showing which root of P a Newton-Raphson iteration will converge to when starting at Zc^\ for the example configuration given in 
Fig. 5. The five positive-sense roots of P are shown by small black triangles, and the unique negative-sense root is shown by a small black circle. 
Each point is colored according to the argument of the root that it converges to, with the central root being shown in white. Darker shades show 
regions that require more iterations to converge. 


gives the iterative equation 


Jn+l) _ J”) I 

+ ip,|2-ip=r 


(63) 


As expected, if P| = 0, then this equation reduces to the stan¬ 
dard complex Newton-Raphson iteration. 

Figure 6 shows a plot of which root the Newton-Raphson 
iteration will converge to as a function of the starting guess 
As is typical for Newton-Raphson iterations of complex 
functions, the plot has a fractal structure, with large basins of 
attraction surrounding each root. However, the plot has some 
distinctly different features to usual Newton fractals [65, 66] 
arising from the vector generalization of the iteration to non- 
analytic functions. In particular, the denominator \P^\'^ — \Pi^ 
featuring in Eq. (63) is zero on a one-dimensional loop of 
points surrounding the central root. Any starting guess that 
approaches this loop will therefore undergo a very large ini¬ 
tial step. In Fig. 6, this loop forms the dividing line between 


the five outer colored basins and the central region. Due to 
the self-similarity of the fractal, the structure surrounding this 
loop is replicated in other parts of the plot. This is in no¬ 
ticeable contrast to the regular Newton fractal for an analytic 
function /(z), where the iteration becomes singular only at a 
zero-dimensional set of points where f'{z) = 0. 


On Figure 6, the object boundary is shown by the dashed 
black line, and it is almost entirely contained within the central 
white region, meaning that a starting guess within the object 
is likely to converge to the collapse point; if the guess is cho¬ 
sen near the center of the object, such as at qo, the iteration 
converges very rapidly and reliably. However the plot also 
indicates that for several small regions inside the object {e.g. 
near the bottom left comer) the Newton-Raphson method may 
converge to one of the exterior roots. 
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1 1.5 


FIG. 7. (a) Time-evolution of a dumbbell-shaped object described by a = \,q\ = —j ^,^3 = —^,5 = and all other Laurent series coefficients 
are zero. The gray curves are plotted at intervals of jL- Positive-sense roots of the function P are shown by triangles, and negative-sense roots 
are shown by circles, (b) The structure of the corresponding function P, with the dashed lines corresponding to contours of |F’(z)| at ” for 
« e N, and the colors corresponding to the argument using the key given in Fig. 5. 


B. Second example: a dumbbell-shaped object dividing in two 

Figure 7(a) shows the dissolution process for the case 
of a long dumbbell-shaped object, where a = I, q\ = 

B = |, and all other Laurent series coefficients are 
zero. In this case, the thin vertical sliver dissolves away leav¬ 
ing two separated fragments. While the system can be time- 
integrated past this point with M(s) remaining non-singular, 
the contour begins to overlap with itself, thus losing physical 
validity. Mathematically, this scenario corresponds to when the 
function g from the unit disk to the physical domain becomes 
multivalued. This is a global property of the function and is 
therefore a different type of finite-time singularity from the 
cusp considered in the previous example. 

The function P for this example is 

P{Zc) =Zc- Btc - qszl + 3qiqiZc - qiZc, 

where tc = I — \q\'^ — 2>\q-i\. The stmcture of P{zc) and its roots 
are plotted in Fig. 7(b). The function P has two negative-sense 
roots in either end of the dumbbell, four exterior positive-sense 
roots, and one positive-sense root on the vertical sliver. This 
example also highlights that the non-analyticity of P signifi¬ 
cantly increases its complexity. The last four terms of P form 
an analytic cubic function in Zc, which could have at most 


three distinct roots, but adding the anti-analytic Zc increases 
the number of roots to seven. 

While more complicated than the previous example, the 
positions of the roots appear to be physically reasonable, with 
one negative-sense root appearing in each end of the dumbbell. 
The central positive-sense root appears to be associated with 
the position where the vertical sliver dissolves. However, close 
inspection reveals that its position is not perfectly aligned 
with the point where the two sides of the object first come 
into contact. Instead, it appears to mark the center of the 
inverted section of the contour at f = L- Figure 8 shows a 
plot of which root the Newton-Raphson iteration converges to, 
depending on the starting guess. For starting points within 
the object, most will converge to the two negative-sense roots 
or the central positive-sense root. The denominator — \Pz!^ 
in Eq. (63) vanishes on two approximate ellipses surrounding 
each negative-sense root. 

While it is not physically valid to simulate the dissolution of 
the object to collapse, this example highlights that the structure 
of P and the position of its roots may be more complicated 
than in the previous example considered, and thus any further 
mathematical analyses would have to take into account this 
possibility. 
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FIG. 8. Plot showing which root of P a Newton-Raphson iteration will 
converge to when starting at Zc^\ for the dumbbell-shaped example. 
The five positive-sense roots of P are shown by black triangles, and 
the two negative-sense roots are shown by black circles. Each point is 
colored according to the argument of the root that it converges to, with 
the central root being shown in white. Darker shades show regions 
that require more iterations to converge. 


C. Third example: transitions in behavior as flow strength is 
altered 

The final example is a three-pronged object given by 
the initial non-zero Lanrent coefficients a = 1, q 2 = 
qs = = ~ 4 (ji ?ii = ~ 1000 ’ ?i4 = Unlike 

the previous two examples, the collapse point equation is diffi¬ 
cult to determine manually due to the large number of Laurent 
series terms that must be considered. However, a computer 
code was written that found it to be 

P{Zc) =Zc-Btc-quzl‘^ +z\'^{-qn + l^uq 2 ) 

+ zl{-q% + llqnq2 + '^^m5-^'iq\Aq2) 

+ 2c (-^5 + 8^8?2 + 11^11 q'5 + 14^14^8 

- 33quq2 - 84^14^2^5 +98^l4g'2) 

+ zli-q2 + 5q5q2 + ^&q5 + llqnq?, 

+ 14^i4g'ii - I2q^ql - 32^iig'2q'5 
-42^14^2^8 -21^140'^ + 22^1 ig'2 
-f 84^i4g'2q'5-35^i4?2)- (64) 

The left panel of Fig. 9 shows the strncture of the solntion 
polynomial when B = 0. Each of the three prongs is surrounded 
by five positive-sense roots, and there is a single negative-sense 
root at the origin. The magnitude of P within the object is small. 


so that most of the object lies within the region |F’(zc)| < 
meaning that an alteration of the flow strength conld alter the 
fnnction’s roots. The right panel shows the function P when 
B = corresponding to a strong flow from the right. In 
this case a new pair of positive-sense and negative-sense roots 
appear on the real axis, resulting in a similar root arrangement 
to Fig. 7(b). 

Figure 10 shows the dissolution process for three different 
cases of B, for a zoomed-in region centered on one of the 
prongs. The top panel shows the case when B = 0, where 
the dissolntion process proceeds normally and the object col¬ 
lapses at the single negative-sense root at the origin. Since the 
physical model given in Eqs. (1), (2), & (3) tends to rapidly 
dissolve sharply curved boundaries, the three prongs of the 
object dissolve rapidly enough that they remain connected to 
each other. 

For the case of Z? = ^ shown in the bottom panel of Fig. 10 
the situation is different. The incorporation of flow into the 
evolution equation of Eq. (11) causes the thin part of the prong 
to dissolve more rapidly than its end, meaning that in this 
case, the object becomes disconnected into two regions. The 
behavior is similar to the previous example, where the object 
boundary overlaps with itself. Att — tc, the object boundary 
loops aronnd the two negative-sense roots and the positive- 
sense root in the same manner as Fig. 7(a). This example 
highlights that only altering flow strength is sufficient to cause 
a transition in the behavior of the dissolution process. 

The transition in behavior is linked with the formation of the 
new roots in Fig. 9 as B is changed from 0 to — However, 
the middle panel of Fig. 10 for an intermediate flow strength of 
B = —^ shows that this transition is more complicated. In this 
case, there is only a single negative-sense root in P. However, 
during time-evolution, the object boundary first overlaps with 
itself, and then the left loop shrinks to zero size, leading to a 
singular solution with an inverted cusp at time t — 0.1550602 < 

tc- 

We carried out a systematic sweep over the flow strengths 
over the range from B = 0toB = —1: initially the object 
collapses to a single point, at B « —0.233 an inverted cusp 
forms, and at B « —0.794 a second negative-sense root forms, 
when the left loop is large enongh to persist nntil f^. This 
result highlights that dissolution process can transition between 
at least three distinct behaviors. Furthermore, the result for 
B = — ^ shows that even if P only has a single negative-sense 
root, the dissolntion process may not be straightforward, and 
may lead to an overlapping boundary or a singular solution. 

Figure 11 shows which roots the generalized Newton- 
Raphson iteration will converge to, for the case of B = 0. The 
plot has an intricate structure and there are many small, distinct 
regions that converge to the central root. The denominator 
— |B||^ in Eq. (63) vanishes on a small loop surronnding 
the central root, and starting guesses near this loop are colored 
in darker shades, indicating that the root-finding algorithm 
takes many iterations to converge. The plot highlights the dif¬ 
ficulty of finding the particular roots of interest in a general 
case. 
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FIG. 9. Plot of the collapse point equation P{zc) for the three-pronged object for (a) S = 0, and (b) B = — ^. The colors correspond to the 
argument of P{z) using the same scale as Fig. 5. The positive-sense roots are shown with black triangles and the negative-sense roots are shown 
with black circles. The thick black line shows the boundary of the object. The thin dashed gray lines are contours of |P(z)| at ^ [n^ + 1) for 
n e Nq. 





VII. CONCLUSION 

In this paper, we studied a model of object dissolution within 
a two-dimensional potential flow, and we created a numerical 
implementation of it that allowed us to simulate the dissolution 
process for arbitrary objects described in terms of a Laurent 
series. The simulations revealed an exact relationship where 
the collapse point Zc is the root of a non-analytic function given 
in the terms of the Laurent coefficients and the flow strength. 
This relationship was subsequently derived analytically, but 
it is unlikely that it would have been discovered without the 
numerical results as a guide. These simulations made use 
of a high-order numerical method, and while these methods 
are often difficult or too computationally expensive to apply 
to real engineering problems, this work demonstrates their 
power in mathematical analysis: the numerical results for the 
collapse point are accurate enough to infer the underlying 
exact relationship with reasonable confidence. There are other 
examples where high-accuracy numerical methods have been 
used for similar purposes, such as demonstrating the existence 
of special solutions to equations [67, 68] or to discovering 
universal behavior [69]. 

The examples of Section VI create some interesting theo¬ 
retical questions for future investigation. We expect that the 
first example of the pentagonal shape (Subsec. VIA) repre¬ 
sents typical behavior for a broad class of objects, where the 
dissolution process is well-defined, the object collapses to a 
single point in finite time, and the collapse point function P(zc) 
has a single negative-sense root. More specifically, we expect 


this to be true for a large class of cases where the qj are small 
in comparison to a, and hence the collapse point function in 
Eq. (58) will be dominated by the Zc term and thus likely to 
have a single anti-analytic root close to the origin. However, 
the second example shows that not all cases may lead to this 
typical behavior, and object may dissolve into multiple frag¬ 
ments, with P{zc) gaining additional negative-sense roots. The 
third example adds a further complication, showing that only a 
minor alteration of the flow strength B can lead to cases where 
dissolution is not well-defined, even though P{zc) still has a 
single negative-sense root. The roots of P{zc) are connected 
to the formation of the finite-time singularities, and give both 
an indication of the formation of local cusps, and potential 
global topological changes. The cusp formation is similar to 
continuous Laplacian growth [23, 29] although a key differ¬ 
ence here is the breakage of symmetry due to the flow. This 
is particularly well illustrated by the third example, where the 
addition of flow breaks the three-fold symmetry and causes 
several transitions in behavior. 

The collapse point results motivate two further questions: 
(A) what the conditions on the initial modes for P to have a 
single anti-analytic root, and (B) what are the conditions on 
the initial modes for the dissolution process to be well-defined 
and for the object to collapse to a single point? If questions A 
and B can be answered, then a further direction would be to 
identify a procedure capable of determining the collapse point 
with absolute certainty. The generalized Newton-Raphson 
method that was introduced in Subsec. VIA is very efficient at 
identifying roots, but it is difficult to determine a priori which 
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FIG. 10. Zoomed-in plot of the dissolution process for the three¬ 
pronged object, for three different values of the flow strength B, 
showing snapshots of the object boundary at intervals of ^ where tc = 
0.1608885. The positive-sense roots are shown with black triangles 
and the negative-sense roots are shown with black circles. For B = 0, 

•JO 

the object collapses to the origin at f = 4- For ^ the object 

boundary overlaps and then a singular solution with a cusp forms 
at t = 0.1550602. For B = —^ the boundary forms an overlapping 
curve at t = fc- 


root it will converge to, and plots of the convergence as a func¬ 
tion of the starting guess exhibit a fractal structure as is typical 
for complex Newton-Raphson iterations. Furthermore, the 
non-analyticity of P(zc) creates some difficulties whereby the 
total number of roots exhibits fundamentally different behavior 
than for analytic functions. An analytic cubic polynomial in 
Zc has exactly three roots (when counted with multiplicity) 
but the addition of a non-analytic Zc as in the second example 
(Subsec. VIB) leads to seven roots, five of which are positive- 
sense and two of which are negative-sense. By using a winding 
argument, considering the curve P{Re‘^) as R ^ we obtain 
— n_ —N where N is the maximum non-zero mode, and n± 
are the number of positive-sense and negative-sense roots. We 
also consider the Newton-Raphson fractals to be interesting 
in their own right, since they have a fundamentally different 
stmcture than typical Newton-Raphson fractals due to the one¬ 
dimensional set of points where the denominator in Eq. (63) 
vanishes. There is an interesting correspondence whereby each 
object has an associated fractal. 

A variety of generalizations to the dissolution model could 
also be explored. The simple form of the right hand side of 
Eq. (11) was based on asymptotic considerations of the concen- 
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FIG. 11. Plot showing which root of the solution polynomial P a 
Newton-Raphson iteration will converge to when starting at Zc^\ for 
the three-pronged object when S = 0. The fifteen positive-sense roots 
of P are shown by black triangles, and the single negative-sense root at 
the origin is shown by a black circle. Each point is colored according 
to the argument of the root that it converges to, with the central root 
being shown in white. Darker shades show regions that require more 
iterations to converge. 

tration profile in the low Peclet number limit, but the numerical 
method could be extended to more complex growth laws where 
higher powers of cos 9 and sin 9 are present. Since the deriva¬ 
tion of the collapse point function is not highly dependent on 
the simple form of Eq. (11), it may be possible to generalize 
this to more complex growth laws as well. Another extension 
is to the case of regular polyhedral objects, which could be 
approximated using many terms in a Laurent series. 

The second and third examples of Subsecs. VIB and VIC 
show that in some cases an object may dissolve into several 
components. In the current numerical method the dissolution 
process cannot be accurately simulated beyond the point where 
multiple fragments form, but it may be possible to extend the 
simulation to this case by using recent advances in conformal 
mapping for multiply connected domains [70-73]. The disso¬ 
lution model is a particularly interesting example, since the 
physical process involves a single domain smoothly transition¬ 
ing into two. We aim to investigate all of these interesting 
directions in future work. 

Einally, we mention the discrete, stochastic analog of this 
problem. In the absence of advection, the diffusion-limited 
erosion of a surface leads to smooth, stable evolution that 
resembles the continuum limit of diffusion-limited dissolu¬ 
tion [21, 35, 42, 43], but to our knowledge this model has 
never been analyzed (or even simulated) to the point where 
the last particle is removed. In the final stages of collapse. 
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discreteness must again become important. Of course, the 
same applies to advection-diffusion-limited dissolution, or 
any other conformally invariant dissolution model [38, 44]. 
We thus leave the reader with an open question: what is the 
probability that a given particle is the last to be removed by 
advection-diffusion-limited erosion of a hnite cluster in a fluid 
flow? 
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Appendix A: Component form of the time-evolution equation 

The numerical method introduced in Section III is based 
upon equating the different sine and cosine components of 
Eq. (14), which is 


( 


N 

\+Ba cos 0 = Re 

ae 

- Y n{bn-icn)e‘"^ 

V 


n=0 


de'^ + Y,ibn + icn)e 

n=0 


Multiplying out these power series yields 


Taking the real component of the bracketed term yields 


N N 
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Collecting terms with factors of sine and cosine yields 
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Equating the terms with different factors of sine and cosine 
yields Eqs. (15), (16), (17), and (18), which together form the 
linear system that is used in the numerical integration method. 
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